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ABSTRACT 

In  this  paper,  we  present  a  new  numerical  formulation  of  solving  the  Boundary 
Integral  Equations  reformulated  from  the  Helmholtz  equation.  The  boundaries  of  the 
problems  are  assumed  to  be  smooth  closed  contours.  The  solution  on  the  boundciry  is 
treated  as  a  periodic  function,  which  is  in  turn  approximated  by  a  truncated  Fourier 
series.  A  Fourier  collocation  method  is  followed  in  which  the  boundary  integral  equation 
is  transformed  into  a  system  of  algebraic  equations.  It  is  shown  that  in  order  to  achieve 
spectral  accuracy  for  the  numerical  formulation,  the  non-smoothness  of  the  integral 
kernels,  associated  with  the  Helmholtz  equation,  must  be  carefully  removed.  The 
emphasis  of  the  paper  is  on  investigating  the  essential  elements  of  removing  the  non¬ 
smoothness  of  the  integral  kernels  in  the  spectral  implementation.  The  present  method 
is  robust  for  a  general  boundary  contour.  Aspects  of  efficient  implementation  of  the 
method  using  FFT  axe  also  discussed.  A  numerical  example  of  wave  scattering  is  given 
in  which  the  exponential  accuracy  of  the  present  numericeil  method  is  demonstrated. 
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tract  No.  NASl-19480  while  the  author  was  in  residence  at  the  Institute  for  t^omputer  Application.s  in 
Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  2.3681. 
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1.  INTRODUCTION 


The  Boundary  Integral  Equation  Method  is  a  powerful  tool  for  solving  certain 
boundary  value  problems.  It  is  particularly  attractive  in  developing  numerical  methods 
since,  when  applied,  it  reduces  the  dimension  of  the  problem  and  often  transforms  a 
problem  in  an  infinite  domain  to  integrals  on  the  finite  boundary  in  which  the  far  field 
radiation  condition  is  satisfied  automatically. 

Numerical  methods  for  the  Boundary  Integral  Equations  have  been  developed 
predominantly  based  on  the  Boundary  Element  Method  (BEM).  In  this  method  the 
physical  bovmdary  is  divided  into  finite  elements  and  integrations  over  each  element  are 
approximated  by  numerical  quadratures.  In  this  way,  the  integral  equation  is  converted 
into  a  system  of  algebraic  equations.  BEM  has  gone  through  a  rapid  advancement  in 
recent  years.  Its  applications  to  the  Helmholtz  equation  are  discussed  in  reference  [2] 
and  the  references  cited  therein.  Other  formulations  of  solving  the  Boundary  Integral 
Equations  in  wave  scattering  and  propagation  have  also  been  proposed  in  the  past, 
including,  for  example,  the  T-matrix  Method[3,4,5]. 

In  this  paper,  we  present  a  spectraJ  method  formulation  for  solving  the  boundary 
integral  equations  reformulated  from  the  Helmholtz  equation.  For  the  Fourier  approx¬ 
imations  used  in  this  paper,  we  assume  that  the  boundary  of  the  problem  is  a  smooth 
closed  contour.  The  fimctions  appearing  in  the  2-D  boxmdary  integral  equation  will 
be  treated  as  periodic  functions,  which  are  in  turn  approximated  with  high  accuracy 
using  truncated  Fourier  series.  The  boundary  integral  equation  is  then  transformed 
into  a  system  of  linear  algebraic  equations.  The  spectral  methods  have  been  known  to 
have  extremely  fast  convergence  rates,  faster  than  any  finite  power  of  1/iV  when  the 
solution  is  infinitely  smooth[6,7]  (where  N  is  the  number  of  collocation  points).  The 
present  numerical  formulation  will  be  seen  to  have  spectral  accuracy. 

It  is  known  that,  although  any  periodic  function  can  be  approximated  by  a  trun¬ 
cated  Fourier  series,  the  rate  of  convergence  depends  on  its  smoothness.  Unfortunately, 
the  integral  kernels  for  the  Helmholtz  equation  axe  not  smooth.  In  particular,  we  note 
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that  the  2-D  Green’s  function,  appearing  in  the  integral  equations,  possesses  a  loga¬ 
rithmic  singiilarity.  Furthermore,  the  normal  derivative  of  the  Green’s  function  also 
contains  a  term  involving  the  logarithmic  function.  It  will  be  seen  that  it  is  critical  to 
remove  the  non-smoothness  of  the  integral  kernel  in  order  to  achieve  fast  convergence 
in  the  spectral  formulation.  In  this  paper  the  non-smoothness  of  the  integral  kernels 
is  subtracted  out  by  using  a  logarithmic  periodic  function  whose  Fourier  expansion  is 
known. 

Foiirier  approximation  methods  for  the  boundary  integral  equations  in  wave  scat¬ 
tering  have  been  proposed  in  the  past  for  simple  geometries.  A  “fast  numerical  method” 
has  been  formulated  by  Bojarski[8]  for  wave  scattering  by  a  circular  cylinder.  It  was 
pointed  out  that  the  bovmdary  integral  equation  for  scattering  by  a  hard  circular  cylin¬ 
der  can  be  solved  easily  and  eflSciently  in  the  Fourier  spectrum  domain  of  the  solution. 
Due  to  the  special  geometry  considered,  a  simple  relation  for  the  Fourier  coefficients 
of  the  solution  and  those  of  the  botmdary  condition  was  derived.  Recently  a  similar 
nmnerical  approach  has  been  applied  by  Schuster[9]  for  a  wave  transmission  problem 
of  concentric  cylinders.  However,  in  these  works  the  singularities  and  non-smoothness 
of  the  integral  kernels  were  not  removed.  Consequently  the  convergence  rates  of  these 
methods  were  not  exponential. 

In  section  2,  a  formulation  of  the  Boundary  Integred  Equation  for  the  Helmholtz 
equation  is  given,  followed  by  a  discretization  using  a  Fourier  collocation  method. 
Essential  elements  of  achieving  spectral  2M:curEM;y  are  investigated.  In  section  3,  a 
numerical  example  is  given  in  which  the  spectrsd  rate  of  convergence  is  demonstrated. 
Section  4  contains  the  conclusions. 


2.  FORMULATIONS 

2.1.  Boundary  Integral  Equations 

Consider  wave  propagation  in  an  interior  or  exterior  domain  with  a  smooth  closed 
boundary  F,  Figure  1.  The  wave  equation  for  a  scalar  function  4>  with  assumed  time 
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dependency  of  e  “*'*  is  reduced  to  the  Helmholtz  equation 


V2^  +  kV  =  0  (1) 

92  ^ 

where  k  is  the  wave  number  and  is  the  2-D  Laplace  operator,  +  -tt-x- 

ox^  d\r 

The  boundary  condition  considered  in  this  paper  is  one  of  the  following  types  : 


Dirichlet :  =  a(f)  On  F 

Neumann  :  C)”  F 

an 

Robin( Impedance)  :  a  ;5— (0  ==  I' 

UTI 

The  Helmholtz  equation  (1)  can  be  reformulated  into  a  Boundary  Integral  Equa¬ 
tion  in  various  ways[10,ll].  Our  purpose  in  this  paper  is  to  demonstrate  the  numerical 
method  using  spectral  approximations  and  the  essential  elements  of  achieving  spec¬ 
tral  accuracy  of  the  numerical  solution.  A  direct  formulation  of  the  Boundary  Integral 
Equation  employing  the  Green’s  Function  will  be  used  here  which  leads  to  an  integral 
equation  involving  <j)  and  its  normal  derivative  on  the  boundary [2]: 


where  —  denotes  the  normal  derivative,  with  the  direction  of  n  being  outward  from  the 
on 

domain  of  wave  propagation,  and  fp  denotes  any  point  on  the  boundeiry.  The  Green’s 
function  G(f,  r"),  whose  form  will  be  given  later,  satisfies  the  following  equation 


V^G  +  K'^G  =  -6{r~f’) 


(3) 


Once  the  values  of  ^  and 


d<f> 

dn 


on  the  botmdary  are  found,  the  solution  ai  any  point. 


Without  loss  of  generality,  let  us  now  assume  the  boundary  curve  F  be  parame¬ 
terized  as 

fr  =  f{e),  0<6<2n  (4) 

where  ^  is  a  non-dimensional  parameter  for  the  boimdary  contour,  not  necessary  the 
arclength  or  some  angle.  In  addition,  the  curve  is  supposed  to  be  simple,  that  is 


dr 

de 


\\dd)  \de)  ^ 


For  simplicity,  we  further  assume  that  f(9)  is  infinitely  differentiable.  Then  the  bound¬ 
ary  integral  equation  (2)  can  be  written  as 


dr 

dd 


de 


(5) 


dr 

d$ 


de.  For  clarity,  the  dependency  on  9  has  been 


in  which  the  line  element  dT  = 
expressed  explicitly  in  (5). 

For  the  2-D  Helmholtz  equation,  the  out-going  Green’s  function  and  its  normal 
derivative  are[12]  : 

9  /« \ 

(6) 


G(e,e')  =  ^Hi'\KR) 


R  •  n 


(7) 


in  which  we  have  used  the  notations  that  R  =  f{9)  —  i^e')  and  R  =  |fi|.  Here  and 
are  the  first  kind  Hankel  functions  of  order  zero  and  one,  respectively. 


i.i.  Spectral  Approximations 


For  a  closed  botmdary  F,  <f>{d)  and  are  periodic  functions  of  0,  for  0  <  0  < 

di 

2x.  Let  <f>  and  -5—  be  approximated  by  truncated  Fourier  series  as  : 

On 


Nl2-\ 
n=-N/2 


in0 


(8) 


4 


(9) 


^^.(6)=  J: 


tn0 


dn 


n=-N/2 

The  pairticular  form  of  the  truncated  Fourier  series  has  been  taken  in  (8)  and 
(9)  for  the  convenience  of  using  discrete  Fast  Fourier  Transform  (FFT)  programs[7]. 
Substituting  (8)  and  (9)  into  (5),  the  boundary  integral  equation  becomes  the  following 


We  note  that  the  two  integrals  in  (10)  are  actually  the  Fourier  coefficients  of 


0(9,6') 


dr 


dd 


QG 

and 


dr 


dd 


the  Fast  Foxirier  Transforms  with  spectral  accuracy. 


respectively.  Our  aim  is  to  evaluate  the  integrals  by 


It  is  clear  that  both  0(6,6') 


dr 


d6 


and  ye.ff) 


df 

d6 


are  periodic  functions  of  6 


and  6'.  They  are  also  infinitely  differentiable  except  ai  6  =  6'  where  H  =  0  in  (6) 
and  (7).  Although  any  periodic  function  can  be  approximated  in  a  truncated  Fourier 
series,  the  rate  of  convergence  depends  on  the  smoothness  of  the  function.  From  this 
consideration,  we  note  that,  first,  0(6,6')  has  a  logarithmic  singularity  at  i?  =  0,  due 

do 

to  the  Hankel  function  of  order  zero  in  (6).  Second,  although  -^(6,6')  can  be  shown 

on 

to  be  a  finite  function,  it  is  not  infinitely  smooth.  It  is  easy  to  show  that  the  function 
R  has  a  discontinuous  derivative  at  0  =  6'.  In  particular,  for  \6  —  6'\  small,  we  have 


R  =  \f(6)  -  f(6')\  = 

=  \6-6'\ 


dr  1  tPr 

3j-(0-»')  +  55^-(0-«T  + 


dv  1  dl^T 

dd  +  2  + 


(11) 


Consequently  any  term  with  am  odd  power  or  logarithmic  fimction  of  R  will  not  be 
infinitely  smooth  and  has  to  be  treated  in  order  to  achieve  spectral  accuracy  in  the 
Fourier  approximation. 


d6 
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To  study  the  singvdarity  in  we  note  that 


0(6,9')  =  +  iYo{KR)) 


in  which  Jq  and  Vo  are  the  zeroth  order  Bessel  functions  of  the  first  and  second  kind, 
respectively.  Using  the  asjrmptotic  series  for  a  small  argument  [13], 

Jo(kR)  =  1 - —  +  — - 

/  n\  2  ,  /  kR\  .  .  ^'y  T  /  ns  (kR)^ 

Yq(kR)  =  -  In  (  —  )  Jo(kR)  H — -Jo{kR)  +  - 

IT  \  2  /  IT  2ir 

we  have,  for  1^  —  ^1  small, 

0(6,6')  =  "-{MKR)  +  iYo(KR)) 

1  /  kR\ 

=  —  —  In  [  j  Jo(kR)  +  (  smooth  terms  ) 

in  which  “smooth  terms”  represents  a  convergent  power  series  containing  only  the 

even  powers  of  R,  and  Jq  is  the  regular  Bessel  function  of  the  zeroth  order[13],  also  a 

power  series  of  R"^.  To  remove  the  logarithmic  singularity  but  preserve  the  periodicity, 

6  ff 

consider  the  periodic  function  In  2sin( — - — ) 

6  =  6'.  Its  Fotirier  series  can  be  foimd  exactly 


In 


2sm(— — ) 


=  -E 


,  which  has  a  logarithmic  singularity  at 

12]  : 

cos  n(6  —  6') 

n=l 


n 


(12) 


To  subtract  out  the  singxilarity  in  0(6,6'),  let 


eCdy)  =  iin 


2sm(— — ) 


Jo(kR) 


Then  the  Green’s  function  is  foimd  as 


0(6,0')  =  0(9,6')  - 


2sm(— 


Jo(kR) 


(13) 


Using  (11),  it  is  easy  to  show  that  0(6,0)  is  periodic  but  finite  for  all  values  of 
the  arguments.  Furthermore  both  0(6,0)  and  Jq(kR)  are  now  infinitely  differentiable 
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functions.  Thus  the  Fotirier  coefficients  of  G{B,6')  —  can  be  computed  according 

du 

to  (13)  with  spectral  accuracy.  For  the  term  involving  the  logarithmic  function  (12), 
convolution  sums  will  be  used. 

We  now  study  the  singularities  in  the  normal  derivative  of  the  the  Green’s  function 

dG 

The  asymptotic  expression  of  H\^\kR)  for  kR  small  hjts  the  following  form[13] 

Ofi 

H[^\kR)  =  Ji{kR)  +  iYi{KR) 

=  Ji{kR)  +  i  ^ - ^  +  —  In  J\{kR)  +  odd  powers  of  kR^ 

2i  2i ,  / kR\  ^  _ 

= - —  -I - In  I  )  Ji{kR)  +  odd  powers  of  kR 

itkR  tt  \  2  / 

in  which  we  have  used  the  fact  that  the  Bessel  function  of  order  unity,  Ji,  is  a  power 
series  of  odd  power  of  kR, 

,  ,  kR  (kRY 

=  +  ••• 

In  addition,  it  can  be  shown  that  R  •  n  =  0[{d  —  6'Y).  Thus  it  follows  that 

)  =  -J  +  -In  (— j  +  odd  powers  of  — 

K  .  ( kR\  ,  ,  .R  •  n  .  ,  . 

=  —  In  I  1  (  smooth  terms  ) 

dG 

It  is  seen  that,  although  is  a  finite  function,  it  does  not  have  smooth  derivatives  due 

on 

to  the  logarithmic  function  appeared  in  the  first  term  shown  above.  For  this  reason, 
its  Fourier  approximation  will  converge  only  at  a  finite  rate  of  1/N^.  To  smooth  out 
the  function  for  computations  by  FFT,  let 


.  ,0-9' 


H(9,0>)  =  -jHY\KRy-^-^ln  2smC-^)  MkR) 


Rn 


Then 


e-e’ 


—(0,e>)=Hi0,9’)  +  —\n\2smC-^)\  MkR) 


R-n 
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R  *  a 

The  functions  and  Ji{kR)—:^  in  (14)  are  now  periodic  and  infinitely 

R 

differentiable.  The  Fourier  coefficients  of  will  be  computed  according  to 

on  du 

(14). 

Following  the  above  analysis,  we  then  express  the  boundary  integral  equation  (10) 
as  follows  : 


N/2-1  /  f2ir  j- 

I  E  E  ^»(/  Te 

_ afirt  _ KT  in  \*'U 


N/2-1 


n=-N/2 


n=-N/2 


K,  ine,  o  ■  /  ^  7  /  Dx  ^ 


R-n  \df\ 


N/2-\  ,  .2n  j=f  7  f 

-  E/»(/ 


n=-N/2 


2ir  fl  —  fl'  dr  \ 

e’^^ln  2sin( — ^ — )  •^o(«^)  ^ 


Now  all  the  integrals  in  (15)  are  in  a  form  which  can  be  evaluated  numerically 
with  spectral  accuracy. 


2.S.  Discretization 

In  what  follows,  a  spectral  collocation  approach  will  be  taken  in  deriving  the 
algebraic  system  for  the  botmdary  integral  equation  (15).  Let  us  introduce  two  sets  of 
discretization  points  (Figure  1) 

ej==2wjlN,  i  =  0,l,...,iV-l 

0;.,  =  2n{j'  +  i)/iV,  y  =  0, 1, ...,  N  -  1 

The  reason  for  using  two  sets  of  discretization  points,  as  will  be  clear  later,  is  that  it 
avoids  the  direct  evaluations  of  G{d,9’)  and  H{6,0')  at  points  where  6  =  6'.  Although 
both  functions  are  finite  there,  their  limits  are  geometry  dependent.  The  current  dis¬ 
cretization  is  robust. 
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We  then  require  that  the  boundary  integral  equation  (15)  be  satisfied  at  collocation 
points  ff  y  =  0, 1, iV  —  1. 

For  convenience  of  discussion,  we  denote  the  following  Fourier  series  approxima¬ 
tions  at  6'j,  for  y  =  0, 1, AT  —  1, 

iV/2-l 

n=-N/2 
N/2-1 

=  E  Vne-'"''  (166) 

n=-N/2 
N/2-1 

=  Pi'ne"”*®  (16c) 

n=-N/2 

5  n/2-1 

tf',  n=-N/2 

The  coefficients  of  the  expansions  are  calculated  by  FFT  (backward  in  the  usual 
sense)  as  follows: 

Si'-  =  jf  ll  %(.»i)  (16“') 

i=0 

’‘i-  =  V  E  §(«y)  (16^) 

j=Q 

Pi'-=jfT,‘‘"’‘lM'cR)].  ^(»i)  (160-) 

>=0  * 

Si•-  =  7^  ^i‘‘"‘‘  ■’,{’‘R)^  §(«i)  (le-?) 

i=o  \e'., 

i' 

In  addition,  we  denote  (12)  as 

In  2sin(^  --  )  =  ^  (16e) 

n=— oo 

where  ao  =  0  and  a„  =  —  ^  for  n  ^  0. 

2  71 1 
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It  is  immediately  seen  that  the  two  integrals  involving  G  and  H  in  (15)  equal  to 
2‘>rgj>n  and  2irhjin  respectively.  The  other  two  integrals  are  obtainable  by  convolution 
stuns.  Specifically,  using  the  definitions  given  in  (16),  we  have 

JV/2-1 

'  (17a) 

m=-N/2 

N/2-1 

d9  =  52  9j'man-me‘^"~'"^76) 

m=-iV/2 

By  substituting  (16)  and  (17)  into  (15),  the  resulting  eJgebraic  system  is  : 

^  N/2-1  N/2-1  N/2-1 

2  52  +  52  (27rh_,v„  +  KUj>n)  =  5Z  -  Vj>n)  (18) 

n=— N/2  n=  — N/2  r»=  — N/2 

for  j'  =  0, 1,2,  ...,N  —  1. 

Equation  (18)  is  easily  cast  into  a  matrix  form 

=  b4!  (19) 

where  4  and  4  are  the  vectors  containing  and  ipn,  respectively.  For  Dirichlet 

A  A 

boimdary  condition,  ^  is  solved  from  (19)  with  $  obtained  firom  the  boimdary  condition 
and  vice  versa  for  the  Neumann  boundary  conditions.  The  Robin  type  impedance 
boundary  condition  can  be  treated  similarly. 

2.4  Evaluation  of  convolution  sums 

The  convolution  sums  in  (17)  require  0{N)  multiplications  for  each  Uj/„  and  Vj'n- 
Thus  the  total  operations  of  obtaining  the  convolution  sums  are  of  order  0{N^  ).  This 
cost  can  be  reduced  considerably  to  0{N'^  logj  N)  by  the  use  of  a  pseudospectreJ  trans¬ 
formation  method  with  de-aliasing  techniques[6,7].  For  completeness,  the  evaluation 
of  (17)  with  a  “padding”  de-aliasing  technique  is  given  below. 

Let  M  >  ZN  and 

ij  =  ,  i  =  0, 1, 2,  ...M  -  1 


Ujv„  =  ^  j  e‘"*ln  2sin( — ^)  Jo(kR)  ^  dd  = 

=  s  i,  '  2“”(-2-)  M’‘R)-g-  Ye 
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Compute  the  following  using  FFT  for  j  =  0, 1, 2,  ...M  -  1  : 


M/2-1 

Aj>j=  ^  amc”"* 

m=  —  M/2 
M/2-1 

pj'i=  E  w--'™*'''"™'' 

m=— M/2 


where 


-  /“ 
Ctm  =  < 

{Vi' 


„  -N  <m<N -I 

0  other 

-iNr/2  <m<N/2-l 

other 


and  form  the  product 


Uj>j  =  Aj>jPj>j 

Then  the  convolution  sum  Ujin  is  the  (backward)  FFT  oiUjij  : 

-  M-l 


i=o 


for  ~N/2  <n<  N/2  —  1.  The  convolution  sum  vj>n  can  be  obtained  identically. 


3.  NUMERICAL  EXAMPLE  :  SCATTERING  OF  A  PLANE  WAVE 
BY  AN  ELLIPTIC  CYLINDER 

The  numerical  method  described  in  the  previous  section  will  now  be  applied  to 
the  problem  of  a  plane  wave  scattering  by  an  elliptic  cylinder.  Exact  solutions  of  wave 
scattering  by  an  elliptic  cylinder  can  be  easily  obtained  in  infinite  series[12].  However, 
numerical  evaluation  of  the  exact  infinite  series  is  not  easily  obtainable  due  to  the 
complexity  of  the  Mathieu  functions  involved.  Our  purpose  here  is  to  demonstrate 
the  exponential  rate  of  convergence  of  the  spectral  method  formulated  in  the  previous 
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section.  The  numerical  results  will  also  be  compared  with  some  asymptotic  values 
available  in  the  literature[14]. 

3.1  Spectral  rate  of  convergence 

Let  the  plane  incident  wave  having  an  incident  angle  a  about  the  x  axis  be 

_  g«ic(ico8ot+y8in  a) 


The  scattered  wave,  satisfies  the  Helmholtz  equation  (1).  The  boundeiry  conditions 

considered  here  are  the  Dirichlet  (soft)  type  <j>  =  —<f>i  or  the  Neumann  (hard)  type 
d<f>  _  d<i>i 
dn  dn  * 

A  simple  parameterization  of  the  elliptic  cylinder  is  given  by 

F:  (x,y)  =  (acos6  ,  bsm0)  ,  0<6<2‘!r  (20) 


where  a  and  b  denote  the  major  and  minor  axis  of  the  ellipse.  The  normal  vector  on  F 
is  ft  =  —(6  cos  0,  a  sin  0)/\/b‘^  cos^  0  +  sin^  9. 

Numerical  results  for  the  Dirichlet  boundary  condition  will  be  presented  first.  For 
this  calculation  the  parameters  have  been  taken  to  be  a  =  2,  6  =  1  and  k  =  27r.  In 
table  I,  the  values  of  the  solution  at  three  selected  points  on  the  boimdary  are  given 

dS 

as  the  number  of  collocation  points  increases.  Listed  are  the  values  of  —  at  0  =  0°, 


dn 

90°  and  180°  on  the  boundary  for  the  incident  angle  a  =  0.  Since  no  exact  value  is 
available,  the  numerical  solutions  are  compeu-ed  with  the  results  obtained  for  N  —  256. 
Clearly  exponential  rate  of  convergence  is  observed. 

We  point  out  that  although  the  parameterization  of  the  ellipse  given  in  (20)  is 
a  convenient  one,  it  may  not  the  be  the  best.  For  example,  the  following  form  of 
parameterization,  which  has  more  uniformly  distributed  collocation  points  than  (20), 

F:  (a:,y)=  cos  (^  +  0.05  sin(2^))  ,  6sin  (^  +  0.05sin(2^))^  ,  0  <  0  <  27r  (21) 

yields  the  values  shown  in  Table  II.  A  somewhat  better  accuracy  for  small  N  is  observed. 
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S.2  Far  field  directivity 

Far  field  scattered  intensities,  computed  as  have  also  been  czdculated.  The 
values  on  an  exterior  point  are  obtained  through  the  boimdary  integral  as 


The  Green’s  function  for  points  lying  outside  of  the  boundary  does  not  have  any  sin¬ 
gularity.  Thus  (22)  is  evaluated  using  FFT  directly.  The  directivities  of  the  scattered 
intensity  for  wave  incident  angle  a  =  0®,  30®,  60®  and  90®  are  shown  in  Figure  3.  For 
this  calcvilation,  a  =  2,  6  =  1  and  k  =  2ir. 

Barakat[14]  gives  some  asymptotic  values  of  long  wave  scattering  of  an  elliptic 
cylinder  at  far  field  for  normal  incidence.  A  comparison  of  the  current  munerical 
results  with  the  asymptotic  values  is  presented  in  Tables  III  and  IV.  For  these  cases, 
a  =  coshpo)  b  =  sinhpo  and  «  =  0.2.  Here  po  corresponds  to  a  value  in  the  elliptic 
coordinates.  The  resxilts  of  the  far  field  scattered  intensity  for  po  =  1-0,  0.6  and  0.2 
are  presented  in  Table  III  and  IV  for  the  Dirichlet  and  Neumann  boundary  conditions, 
respectively.  For  the  present  long  wave  scattering,  V  =  64  has  been  used  for  eill  the 
calculations.  The  numerical  computation  agrees  with  the  aisymptotic  estimation. 

Finally  we  note  that  although  the  exterior  scattering  problem  is  uniquely  deter¬ 
mined,  the  direct  formulation  of  the  Boundary  Integral  Equation  used  here  would  not 
yield  a  unique  solution  when  the  wave  number  «  coincides  with  the  eigenvalues  of  a 
correspondent  interior  homogeneous  problem.  On  the  other  hand,  this  numerical  diffi¬ 
culty  is  well  tmderstood  and  remedies  are  readily  available[10,ll].  For  example,  vmique 
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solution  can  be  obtained  in  a  combination  of  single-  and  double-layer  formulation  of  the 
Boimdary  Integral  Equation[15].  The  modifications  of  the  present  spectral  implemen¬ 
tation  to  other  formulations  of  the  Boimdary  Integral  Equation  are  straightforward. 


4.  CONCLUSIONS 

In  this  paper,  a  spectral  method  of  solving  the  Boundary  Integral  Equation  is 
presented.  It  is  shown  that  the  integral  kernels  for  the  Helmholtz  equation  contain 
singular  terms  that  have  to  be  removed  to  achieve  the  spectral  accuracy.  Detailed 
numerical  implementation  of  a  Fourier  collocation  formulation  has  been  given.  The 
non-smoothness  of  the  integration  kernels  is  subtracted  out  by  using  a  logarithmic 
function  whose  Fourier  expansion  is  known.  The  numerical  formulation  presented  here 
preserves  the  spectral  accuracy  and  yields  an  exponential  rate  of  convergence. 

Compared  to  the  Boundary  Element  approaches,  the  Spectral  Boundary  Integral 
Equation  Method  presented  in  this  paper  would  yield  matrices  of  much  smaller  size 
since  the  latter  requires  far  fewer  points  to  achieve  the  same  accuracy.  This,  of  course, 
reduces  the  complexity  and  the  cost  of  solving  the  resultant  algebraic  system.  Since 
both  methods  will  result  in  dense  matrices,  it  appears  that  the  spectral  formulation  is 
more  advantageous  for  problems  with  smooth  geometries. 
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Values  of 


d<i> 

dn 


TABLE  I 

for  the  Dirichlet  boundary  condition  at  selected  points  on  the  boundary. 

The  boundary  is  given  by  (20). 
a  =  2,  6  =  1,  K  =  27r,  incident  angle  a  =  0°. 


N 

II 

O 

o 

6  =  90" 

O 

O 

00 

II 

Relative  Error 

32 

6.178009567 

2.376671291 

7.134225119 

0.04 

48 

6.379285016 

2.175324633 

6.914464114 

0.2  xlO-'* 

64 

6.379334010 

2.175259123 

6.914387795 

0.1x10-'^ 

80 

6.379333962 

2.175259114 

6.914387773 

10"^ 

96 

6.379333961 

2.175259114 

6.914387772 

10"" 

256 

6.379333961 

2.175259114 

6.914387772 

Values  of 


dn 


TABLE  li 

for  the  Dirichlet  boundary  condition  at  selected  points  on  the  boundary. 

The  boundary  is  given  by  (21). 
a  =  2,  6  =  1,  K  =  27r,  incident  angle  a  =  0°. 


N 

II 

o 

o 

0  =  90" 

o 

00 

II 

Relative  Error 

32 

6.402635098 

2.165098187 

6.899520927 

0.005 

48 

6.379355527 

2.175263759 

6.914393819 

0.3x10"® 

64 

6.379334276 

2.175259159 

6.914387893 

0.5x10"'^ 

80 

6.379333966 

2.175259114 

6.914387774 

0.5x10"^ 

96 

6.379333961 

2.175259114 

6.914387772 

10"" 

256 


6.379333961 


2.175259114 


6.914387772 


TABLE  III 

Values  of  r<f>^  at  far  field  for  the  selected  angles.  Dirichlet  type 
soft  wall  botindary  conditions  are  applied. 
a  =  coshpo,  h  —  sinhpo?  «  =  0-2  and  incident  2Lngle  a  =  90°. 


Far  Field  Angle 

Numerical 

Asymptotic[14] 

po  =  1.0,  e  =  0°,180° 

1.769057417 

1.774 

po  =  1.0,  6  =  90° 

2.137413140 

2.113 

po  =  1.0,  d  =  270° 

1.538098511 

1.491 

/io  =  0.6,  6  =  0°,  180° 

1.364415306 

1.367 

Hq  =  0.6,  e  =  90° 

1.506576730 

1.484 

po  =  0.6,  d  =  270° 

1.287260683 

1.258 

/io  =  0.2,  6  =  0°,  180° 

1.062775602 

1.064 

po  =  0.2,  6  =  90° 

1.107062467 

1.087 

/io  =  0.2,  6  =  270° 

1.062022484 

1.041 
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TABLE  IV 

Values  of  r<f>^  at  far  field  for  the  selected  angles.  Neumann  type 
hard  wall  boundary  conditions  are  applied. 
a  =  cosh^Oi  b  =  sinh/io,  «  =  0.2  and  incident  angle  a  =  90°. 


Far  Field  Angle 

[ 

Numerical 

Asymptotic[14] 

fiQ  =  1.0,  e  =  0",180" 

0.009175011010 

0.00884 

fio  =  1.0,  e  =  90" 

0.02258893859 

0.0216 

tio  =  1.0,9  =  270" 

0.1110188572 

0.1113 

Po  =  0.6,  6  =  0",  180" 

0.001676516587 

0.00165 

^0  =  0.6,  9  =  90" 

0.007209616977 

0.0071 

fio  =  0.6,  9  =  270" 

0.02745548914 

0.0273 

fiQ  =0.2,  fl  =  0",180" 

0.0001288604229 

0.00013 

fiQ  =  0.2,  9  =  90" 

0.003710427872 

0.0037 

po  =  0.2,  9  =  270" 

0.006989684179 

0.0069 
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Figure  1.  Schematic  of  discretization  of  the  boimdary.  6  is  the  parameterizing  variable 
for  the  boimdary  contour,  n*  and  fZ,  indicate  the  direction  of  the  normal  vector  for  an 
exterior  and  interior  problem  respectively. 
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d<f> 

Figure  2.  Solution  -5—  on  the  surface  of  the  ellipse  for  the  Dirichlet  boundary  condition. 
cfn 

a  =  2,  6  =  1  and  k  =  27r.  Plane  wave  incident  angle  a  =  0".  The  parameterization  of 
the  ellipse  is  that  of  (21). 
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Figure  3.  Far  field  directivities  of  the  scattered  wave  for  indicated  incident  angles, 
a  =  2,  6  =  1  and  «  =  27r,  with  Dirichlet  boundeury  condition. 
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